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Abstract 

Developmental conservation among related species is a common generalization known as von Baer's third law and implies that early 
stages of development are the most refractory to change. The "hourglass model" is an alternative view that proposes that middle 
stages are the most constrained during development. To investigate this issue, we undertook a genomic approach and provide 
insights into how natural selection operates on genes expressed during the first 24 h of Drosophila ontogeny in the six species of the 
melanogaster group for which whole genome sequences are available. Having studied the rate of evolution of more than 2,000 
developmental genes, our results showed differential selective pressures at different moments of embryogenesis. In many Drosophila 
species, early zygotic genes evolved slower than maternal genes indicating that mid-embryogenesis is the stage most refractory to 
evolutionary change. Interestingly, positively selected genes were found in all embryonic stages even during the period with the 
highest developmental constraint, emphasizing that positive selection and negative selection are not mutually exclusive as it is often 
mistakenly considered. Among the fastest evolving genes, we identified a network of nucleoporins (Nups) as part of the maternal 
transcriptome. Specifically, the acceleration of Nups was driven by positive selection only in the more recently diverged species. 
Because many A/i/ps are involved in hybrid incompatibilities between species of the Drosophila melanogaster subgroup, our results link 
rapid evolution of early developmental genes with reproductive isolation. In summary, our study revealed that even within functional 
groups of genes evolving under strong negative selection many positively selected genes could be recognized. Understanding these 
exceptions to the broad evolutionary conservation of early expressed developmental genes can shed light into relevant processes 
driving the evolution of species divergence. 
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Introduction 

During ontogeny, most hierarchical features are a conse- 
quence of the timing of developmental events. Indeed, as 
later events depend on earlier ones, developmental constraints 
during embryonic stages are widespread (Carroll et al. 2001). 
As a consequence, genes involved in early developmental pro- 
cesses are expected to be under strong negative selection to 
prevent deleterious cascading effects (Roux and Robinson- 
Rechavi 2008; Artieri et al. 2009). In Drosophila, it has been 
recently shown that embryonic genes evolve at a slower pace 
than postembryonic and adult expressed genes (Artieri et al. 
2009). However, the pattern of early conservation has not 
been supported when embryonic-specific analyses were 



carried out (Davis et al. 2005; Cruickshank and Wade 2008; 
Kalinka et al. 2010). 

Drosophila development is characterized by a fast segmen- 
tation process. Segment determination starts very early in em- 
bryogenesis, when approximately 3h after fertilization the 
position and identity of all body structures are determined 
simultaneously during the blastoderm stage (Foe and Alberts 
1983). Although not all genes involved in segmentation have 
an exclusive timing of expression, two hierarchical regulatory 
layers can be identified: An initial phase of maternal compo- 
nent specification and a succeeding phase involving more 
complex and interactive zygotic gene expression (Schroeder 
et al. 2004). Genes in early layers (maternal genes) regulate 
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the expression of genes in subsequent layers (gap, pair-rule, 
segment polarity, and Hox genes) but not vice versa (Jaeger 
2009). In addition, there is cross-regulation among genes in 
the same hierarchical layer (Manu et al. 2009). Understanding 
these two phases and their connections is necessary to recon- 
struct the evolution of the embryonic system (Wilkins 2002). 

The developmental stage that is most refractory to evolu- 
tionary change is commonly known as the phylotypic stage 
(Raff 1996). Based on genome-wide expression comparisons 
(Kalinka et al. 2010) and sequence evolution analyses (Davis 
et al. 2005; Cruickshank and Wade 2008), the initiation of 
organogenesis during the burst of expression of segment po- 
larity and Hox genes appears to be the Drosophila phylotypic 
stage. Instead, earlier embryonic stages, including the mater- 
nal component of segmentation, have markedly diverged 
within and among species (Galis et al. 2002). Thus, given 
that the highest constraint takes place during middle embryo- 
genesis, a developmental hourglass model likely reflects 
Drosophila embryonic evolution (Raff 1986). 

Despite the strong developmental constraint across phylo- 
typic stages, cases of rapid evolution at embryonic expressed 
genes were identified. Specifically, the rapid evolution 
reported for both maternally and zygotically expressed Hox 
and Hox-derived genes challenges the view of general conser- 
vation of embryonic genes (Barker et al. 2005; Casillas et al. 
2006). Whether these cases of rapid evolution are driven by 
positive selection (PS) or just relaxation of selective constraints 
(RSCs) remains unknown. Moreover, although the incidence 
of fast evolving genes expressed in early development is ex- 
pected to be low, there is a lack of studies searching for the 
signatures of PS and hence fast adaptive change in embryonic 
genes. 

Taking advantage of fine time-course gene expression in- 
formation (Hooper et al. 2007), the recent burst of Drosophila 
species whole genomes sequences (Clark et al. 2007) and the 
development of powerful statistical and bioinformatic tools 
(Yang 2003, 2007), here we investigate the evolution 
of tightly regulated groups of embryonic genes in the 
Drosophila melanogaster species group. Specifically, we stud- 
ied the evolutionary rates of genes involved in the three major 
embryonic groups recognized by Hooper et al. (2007): 
Maternal, early zygotic, and late zygotic genes, all predomi- 
nantly expressed in a stage-specific fashion during embryo- 
genesis. In addition, we performed specific maximum 
likelihood (ML) tests to distinguish true cases of PS from 
likely cases of RSC (Serra et al. 201 1). Importantly, our geno- 
mic-scale study allowed us not only to dissect the rapidly evolv- 
ing fraction of the Drosophila embryonic transcriptome but 
also to contrast the incidence of positively selected genes 
(PSG) at the phylotypic stage with less constrained periods 
of development. After identifying rapid evolving genes ex- 
pressed during embryogenesis, a functional analysis of this 
particular fraction of the genome was performed to shed 
light into the processes involved in developmental adaptation. 



Having studied the evolution of more than 2,000 embry- 
onic genes, our results are in agreement with the hourglass 
model of evolution. However, the incidence of PSGs was 
homogeneous across embryonic development. Among the 
fastest evolving genes, we identified a network of nucleopor- 
ins genes (Nups) as part of the maternal transcriptome. 
However, these rapidly evolving genes exhibited signatures 
of adaptive evolution only in the most recently diverged spe- 
cies of the melanogaster species group studied. Because Nups 
have been shown to be involved in hybrid incompatibilities 
(Tang and Presgraves 2009; Sawamura et al. 2010), we dis- 
cuss the possible role of nuclear pore-related developmental 
processes in species divergence. 

Materials and Methods 

Data Acquisition 

Information on the timing of gene expression was obtained 
from a previous survey of embryonic gene expression that 
performed a time-course genome-wide microarray analysis 
(Hooper et al. 2007). The survey consisted of an extensive 
analysis of the fly transcriptome obtained in 30 time points, 
uncovering the entire 24-h period in which the fertilized egg 
develops into a first-instar larva. By applying convolution 
methods (e.g., common "sharp" transcript changes among 
genes), the authors identified three major categories for which 
all transcript levels increase and/or decrease within a certain 
time interval, suggesting a common mode of regulation. The 
first group includes highly expressed genes that encode ma- 
ternal transcripts that show a subsequent decrease in expres- 
sion by 12 h after egg-laying. The second group consists of 
early zygotic genes with high transcription levels starting at 2- 
3 h after egg laying (embryo stage 5) and that later decrease in 
midembryogenesis. Finally, late zygotic genes encode tran- 
scripts for which we only observe an increase in expression 
starting 1 3 h after egg laying (embryo stage 1 6-1 7), maintain- 
ing high expression levels till the end of embryogenesis. The 
total number of genes for each embryonic stage included in 
our analysis along with the proportion of sex-biased expres- 
sion is given in table 1 . Gene orthology relationships among 
the six Drosophila species studied could be ascertained for only 
62% of the genes studied by Hooper et al. (2007). Thus, 999 
out of 1 ,534 of the maternal genes (class I), 496 out of 792 of 
the early zygotic genes (class II), and 597 out of 1,053 of the 
late zygotic genes (class III) could be included in the present 
report. Classes refer to the nomenclature used by Hooper 
et al. (2007). 

Evolutionary Rate Estimation 

Coding sequences (CDS) data of embryonic genes were ob- 
tained from the genomes of 1 2 Drosophila species available at 
www.flybase.org (last accessed November 18, 2013; Clark 
et al. 2007). CDSs were aligned with MUSCLE (Edgar 2004) 
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Table 1 

Groups of Embryonic Genes and Information of Sex-Biased Expression 



Groups of Genes 


No. of 


Female-Biased 


Male-Biased 


Unbiased 


Unclassified 




Genes 


Genes 


Genes 


Genes 


Genes 


Maternal genes 


999 


805 (80.6%) 


51 (5.1%) 


99 (10%) 


44 (4.3%) 


Early zygotic genes 


496 


151 (30.4%) 


87 (17.5%) 


246 (49.6%) 


12 (2.5%) 


Late zygotic genes 


597 


31 (5.2%) 


258 (43.2%) 


237 (39.7%) 


71 (11.9%) 



using predicted amino acid sequences as templates. Aligned 
columns containing gaps were removed. We only included in 
the analysis the six species of the D. melanogaster group be- 
cause saturation in silent site divergence outside the D. mela- 
nogaster speoes group precludes the use of all 12 genomes 
(Larracuente et al. 2008). For each gene for which orthology 
could be confidently determined, we calculated the ratio of 
nonsynonymous (d/V) to synonymous (d5) substitutions rates 
(co) for each species. Estimates of co were obtained applying a 
free ratio branch ML model using CodeML program of the 
PAML 4 package (Yang 2007). Values of evolutionary rates 
(d5, d/V, and co) for all genes analyzed are shown in supple- 
mentary table S1, Supplementary Material online. 

PS was evaluated using two different branch-site models (A 
and A1) also implemented in CodeML (Yang 2007). Branches 
in the phylogeny were defined a priori as foreground and 
background lineages. Under these models, only the fore- 
ground lineage may contain events of PS. The test was per- 
formed independently for each species by marking its 
corresponding terminal branch, this mark indicates that a spe- 
cific evolutionary model may be applied to the flagged branch 
(either evolving under three main classes of evolutionary rate 
&> 0 < 1 , = 1 , &>2a/b > 1 ; or in the case of A1 model, only 
under co 0 or co^). Because the compared models are nested, 
likelihood ratio tests were performed and likelihood ratio tests 
statistics {2A€ = -2 [In (likelihood for null model) -In (likeli- 
hood for alternative model)]} were posteriorly transformed 
into exact P values using the pchisq function of the R statistical 
package. Likelihood ratio tests were performed using a x 
distribution with df = 2 for Test A and df = 1 for Test A1, 
which have been shown to be conservative under conditions 
of PS (Zhang et al. 2005). P values derived from PS analyses 
were false discovery rate-adjusted using the method of 
Benjamini and Hochberg (1995). In contrast to the statistical 
behavior of previous branch-site tests, the methodology pro- 
posed by Zhang et al. (2005) represents an improvement in 
this kind of test based on the comparison of the ML of differ- 
ent evolutionary models. This approach has proved to be 
able to successfully differentiate PS from RSCs and weak PS. 
For further details of the parameters used in ML models, see 
Lavagnino et al. (2012). All adjusted P values of branch- 
site models are reported in supplementary table S2, 
Supplementary Material online. In addition, supplementary 
table S3, Supplementary Material online contains all genes 
evolving under PS in all species studied. 



Heat Map Construction and Clustering Analysis 

We employed Ward's (1963) method for clustering analy- 
sis using 6N/6S values from free ratio branch ML model. 
Clusters were defined using the 99.9 percentile of the dis- 
tance matrix as threshold value. d/V/d5 values upper than 1 
were considered equal to 1 to avoid bias toward infinite 
values. 

Gene-Set Selection Analysis 

We performed a gene-set enrichment analysis employing the 
program BABELOMICS (Al-Shahrour et al. 2007) to study the 
association of specific Gene Ontology (GO) terms to fast evolv- 
ing genes in each particular stage studied (maternal, early zy- 
gotic, and late zygotic). For this reason, all genes were ranked 
according to the co value in D. melanogaster and looked for 
blocks of functionally related genes in the group of the fast 
evolving genes. The program also corrected P values for mul- 
tiple adjustment effects by false discovery rate. 

Identification of Networks of PSG 

We searched for networks of interacting PSG using R-spider 
(Antonov et al. 2010). This tool determines whether interac- 
tions between input genes are greater than expected by 
chance. Input data were the complete list of genes classified 
as PSG in all Drosophila species studied. Finally, employing 
STRING database (Szklarczyk et al. 2011), the network was 
built with all Nups and interacting PSG allowing one missing 
node. 

Temporal Specificity of Nups across Development 

We adapted the tissue specificity index, r. to investigate how 
narrow Nups network genes expression is across 
development: 

N I c 

El _ logf,- 
logEmax 



where N is the number of stages being compared, E,- is the 
expression in stage /, and E max is the maximum expression 
reached by the gene across stages (Yanai et al. 2005). In our 
case, we employed expression data from time-course 
genome-wide microarray analysis of Graveley et al. (201 1). r 
ranges from 0 to 1, with larger r values indicating greater 
temporal specificity. 
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Results 

Rates of Evolution in Coregulated Groups of Genes 
during Embryogenesis 

ML estimates of 6N/6S of 2,092 embryonic genes (table 1 
and supplementary table S1, Supplementary Material online, 
for detailed information) using a free ratio branch model re- 
vealed significant differences among maternal, early zygotic, 
and late zygotic groups of genes in all Drosophila species, 
except D. sechellia (fig. 1 ). In D. melanogaster and D. simulans, 
maternal genes evolved significantly faster than early 
zygotic genes, while maternal genes evolved faster than 
late zygotic genes in D. melanogaster (fig. 1). By contrast, 
late zygotic genes evolved significantly faster than maternal 
and early zygotic genes in D. erecta, D. yakuba, and D. ana- 
nassae (fig. 1). 

A clustering analysis using d/V/d5 values of each gene from 
the six species identified four different groups of genes (fig. 2). 
The largest cluster included genes with the lowest 6N/6S 
values and clusters 2, 3, and 4 contained genes with fast 
evolutionary rates. Interestingly, maternal, early zygotic, and 
late zygotic genes were not randomly distributed across 
clusters (x 2 =10.9, P= 0.0026). Specifically, slow evolving 
genes (cluster 1) were enriched in early zygotic genes. 
Such pattern indicates a shared signature of purifying selec- 
tion during middle embryogenesis across species. On the 
contrary, fast evolving genes grouped in clusters 3 and 4 ex- 
hibited lineage-specific acceleration in D. sechellia and 
D. simulans, respectively. Finally, a common across-spe- 
cies signature of rapid evolution was detected for genes of 
cluster 2. 

A gene-set selection analysis (Serra et al. 2011) detected 
functional sets of rapidly evolving genes in each one of the 
three embryonic stages of D. melanogaster development 
(table 2). Notably, we found that rapidly evolving genes of 
maternal and early zygotic expression are enriched in compo- 
nents of intracellular membranes such as "pore complex" and 
"organelle envelope" in maternal expressed genes and "in- 
tracellular membrane-bound organelle" in the early zygotic 
genes (table 2). 

It is known that gene expression level and genomic location 
are among the most important factors affecting evolutionary 
rates in Drosophila species (Clark et al. 2007). Thus, we inves- 
tigated the distribution of the genes included in the three 
major embryonic groups among chromosomes and the rela- 
tionship between gene expression level and evolutionary rates 
to rule out the possibility that these factors other than devel- 
opmental timing of expression are shaping the differential 
evolutionary rates observed. First, we found that the genes 
included in the three major embryonic groups are randomly 
distributed in the genome (x 2 =7.4, P= 0.285). Second, 
though a regression analysis of the entire set of genes used 
in our study indicates that highly expressed genes evolved 
more slowly than less expressed ones (F 1#174 o= 18.09, 



P< 0.0001), a Kruskal-Wallis analysis of variance (ANOVA) 
showed that differences in the level of gene expression 
among genes involved in the different embryonic stages 
were not significant (P= 0.085). Thus, we can argue that al- 
though gene expression level and genomic location are factors 
known to influence gene sequence evolution in Drosophila 
species these factors do not affect the comparisons performed 
in this study. 

PSG in the Embryonic Transcriptome 

We searched for cases of PSG employing the tests developed 
by Zhang et al. (2005), which permit to distinguish between 
cases of PS from false positives due to RSCs (or weak signals of 
PS). Interestingly, PSGs were found in all stages and species 
studied (supplementary table S3, Supplementary Material 
online). Surprisingly, a comparison of the incidence of PS be- 
tween embryonic and adult stages in the six species of the 
D. melanogaster group analyzed revealed a higher proportion 
of PSG in embryonic transcriptomes in D. simulans, D. sechel- 
lia, D. erecta, and D. ananassae as shown in table 3. 
Nevertheless, it is worth mentioning that evolutionary rates 
of embryonic genes were significantly lower than male-adult 
expressed genes for all Drosophila species (fig. 3), indicating 
that the general slower pace of evolution is independent of 
the high incidence of PS in embryonic genes of the aforemen- 
tioned species. 

Network of Embryonic PSG 

After the identification of PSG in the embryonic stages of 
the six species of the D. melanogaster species group se- 
quenced so far, we searched for networks of interacting 
PSG using the program R-spider (Antonov et al. 2010). 
Interestingly, only a single network of PSG was identified 
in the maternal transcriptome (P= 0.025). This module con- 
sists of a set of nucleoporins (Nups) genes that encodes 
proteins involved in the structure of the nuclear pore com- 
plex. Within this network, we found cases of PSG in all 
species studied (table 4). To further investigate the evolution 
of the Nups network, we divided the sample of embryonic 
genes into two groups: The Nups network and the rest of 
embryonic genes and compared their rate of evolution in 
each species. We found that the Nups network evolved 
faster than the rest of the embryonic genes in D. sechellia, 
D. simulans, and D. melanogaster but not in D. erecta, D. 
yakuba, and D. ananassae (fig. 4). Moreover, to identify at 
which point in the phylogeny started the acceleration of 
Nups evolution, we compared the nonsynonymous substitu- 
tion rates of Nups in each lineage with the rates calculated 
for the respective most recent common ancestors in the 
Drosophila phylogeny. These tests also revealed a significant 
increase in the pace of the evolution of Nups though only in 
D. sechellia (fig. 4; Wilcoxon-matched pairs test: A/=15, 
Z=2.7, P<0.01). Interestingly, these results are in 
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Fig. 1. — Box plots show the distribution of d/V/cLS values for maternal (M), early zygotic (EZ), and late zygotic (LZ) genes for the six species of the 
D. melanogaster group. Each box extends from the first to the third quartile, with the line in the middle of the box indicating the median. Asterisks show a 
significant difference in d/V/d5 between different developmental stages in a Kruskal-Wallis test. For D. melanogaster. M versus EZ, P= 0.020, M versus LZ, 
P= 0.026; for D. simulans: M versus EZ, P= 0.01 9; for D. yakuba: M versus LZ, P < 0.001 , EZ versus LZ, P= 0.006; for D. erecta: M versus LZ, P < 0.001 , EZ 
versus LZ, P < 0.001 ; for D. ananassae: M versus EZ, P < 0.001 , M versus LZ, P < 0.001 , EZ versus LZ, P= 0.044. To easier the lecture, the plot was truncated 
for d/V/d5 values >0.5. 
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Fig. 2. — Heat map using d/V/d5 values for the six species of Drosophila studied. We employed ward method for clustering analysis using dN/dS values 
derived from free ratio branch ML model. 



agreement with a recent report showing nucleoporins as a 
preferential target of PS in D. mauritiana, in which D. sechel- 
lia is a member of the simulans clade (Nolte et al. 2013). All 
in all, these results suggest that rate of evolution of Nups is 
only accelerated in recently diverged species of the melano- 
gaster group but not in older lineages. 



Temporal Specificity of Nups across Development 

Because a positive association between temporal specificity 
of gene expression and evolutionary rate has been reported 
(Artieri et al. 2009), we analyzed whether such pattern 
occurred for Nups network genes. For this reason, we 
calculated Nups temporal specificity across development 
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Table 2 

Overrepresentation of GOs in Fast Evolving Genes 



Groups of Genes 



Fast Evolving Genes 



Maternal genes 



Immune response (GO:0006955) 
Response to biotic stimulus (GO:0009607) 
Pore complex (GO:0046930); organelle 
envelope (GO:0031967) 



Early zygotic genes 



Intracellular membrane-bound organelle 
(GO:0043227) 



Late zygotic genes 



Cellular protein metabolic process 

(GO:0044237) 
Peptidase activity (GO:0008233) 
Intracellular organelle part (GO:0044446) 



Note. — Red refers to GO biological function terms, blue to GO molecular 
function terms, and green to GO cellular component terms. 



Table 3 

Comparison of the Incidence of PSGs in Embryonic and Adult 
Expressed Genes of the Six Drosophila Species 



Species 



Embryonic PSGs Adult PSGs 



D. melanogaster 


13 (0.62%) 


9 (0.9%) 


ns 


D. simulans 


61 (2.92%) 


10 (1%) 


10 9*** 


D. sechellia 


47 (2.25%) 


6 (0.6%) 


10.8*** 


D. yakuba 


9 (0.43%) 


3 (0.3%) 


ns 


D. erecta 


14 (0.67%) 


1 (0.1%) 


4.5* 


D. ananassae 


62 (2.96%) 


7 (0.7%) 


15.7*** 



Note. — Number of embryonic expressed genes = 2,092; number of adult ex- 
pressed genes = 993. *P<0.05, **P<0.01, ***/>< 0.001. ns: not significant 



using the tissue specificity index (Yanai et al. 2005). The 
results of this analysis revealed that though Nups exhibited 
intermediate temporal specificity values with a mean of 
0.27, they showed the highest expression during embryo- 
genesis (fig. 5). In this sense, even if we cannot rule out the 
possibility of Nups function during nonembryonic stages of 
development, the highest expression level during embryo- 
genesis is a reliable indicator of an important embryonic 
function. 



Discussion 

The present evolutionary genomics study demonstrates that 
coexpressed groups of genes across embryogenesis are sub- 
ject to differential selective pressures. Though we confirm the 
pervasive role of negative selection in early development, we 
identified a large number of PSG as part of the embryonic 
transcriptome. Remarkably, a fast evolving network of nucleo- 
porins stands as an island of rapid embryonic evolution. 
Altogether, our findings highlight that despite being part of 
the stage with strongest developmental constraint across 
Drosophila ontogeny, many embryonic genes show rapid evo- 
lutionary change. 



Hourglass-Type Embryonic Evolution 

Having studied the evolution of 2,092 embryonic genes, we 
may conclude that early zygotic genes, but not maternal and 
late zygotic genes, are subject to the strongest evolutionary 
constraints during embryogenesis in D. melanogaster and 
D. simulans (fig. 1). Moreover, when analyzing embryonic 
genes evolution across species, we found a shared signature 
of purifying selection mainly in early zygotic expressed genes 
(fig. 2). These results do not support the developmental con- 
straint hypothesis and states that genes involved in early de- 
velopmental processes are under strong negative selection to 
prevent deleterious cascading effects (Artieri et al. 2009). On 
the contrary, our results are in agreement with the embryonic 
"hourglass" model that posits the onset of segmentation as 
the Drosophila phylotypic stage (Raff 1986). Such distribution 
of the rate of evolutionary change of genes expressed at dif- 
ferent moments in embryonic development can be mainly 
explained by two independent causes. First, the strongest re- 
fractory period to evolutionary change across ontogeny takes 
place when organogenesis begins during the early zygotic 
stage. In addition, the early zygotic transcriptome exhibits 
the highest number of protein interactions during embryogen- 
esis (Hooper et al. 2007), a fact that may also contribute to the 
strong functional constraint. Second, the maternal trancrip- 
tome is mainly composed of genes with female-biased expres- 
sion (table 1). Sex-biased and reproduction-related genes are 
among the fastest evolving genes in animal genomes (Meisel 
201 1 ; Assis et al. 201 2). Thus, on the one hand, female-biased 
gene expression may drive maternal transcriptome accelera- 
tion, and, on the other hand, the relevance of developmental 
processes during the onset of segmentation imposes a re- 
striction on evolutionary change in early zygotic genes. 
Interestingly, another hourglass pattern was recently reported 
in other stages of Drosophila development. In effect, greater 
conservation of gene expression levels during the pupal stage 
was found in comparison with third-instar larvae and adult 
stages in species of the D. melanogaster subgroup and inter- 
specific hybrids (Artieri and Singh 2010). Such hourglass pat- 
terns are not unexpected because many genes share a 
biphasic expression pattern during development in the early 
embryo and later in the pupal stage. Strikingly, rounds of ex- 
tensive organogenesis with regulatory conservation are shared 
by embryogenesis and the pupal stage (Arbeitman et al. 
2002). We may add that though at the level of protein- 
coding sequence and expression patterns (Kalinka et al. 
2010) middle embryogenesis seems to be the most con- 
strained stage of development we cannot assume that at 
level of regulatory sequences an hourglass model may also 
fit as well. On the contrary, in vertebrates, the evidence 
points to the validity of the hourglass model only for regula- 
tory regions but not for protein-coding genes (Piasecka et al. 
2013). 
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Fig. 3. — Box plots show the distribution of d/V/d5 values for embryonic and adult expressed genes for the six species of the D. melanogaster group. Adult 
expressed genes evolved significantly faster than embryonic genes in all species (Kolmogorov-Smirnov test, P< 0.001). Adult expressed genes were 
represented by clusters 18-20 from Graveley et al. (201 1). To easier the lecture, the plot was truncated for d/V/d5 values >0.8. 



Table 4 

PS at Nucleoporins and Interacting Partners in Drosophila Species 



Gene 


Interactions 


Positive Selection 


Pen 


9 


D. simulans and D. yakuba 


CG4887 


2 


D. simulans 


CG7185 


3 


D. simulans 


Nup154 


19 


D. simulans 


Nup160 


13 


D. simulans 


Cpsf160 


15 


D. simulans 


dgt5 


3 


D. sechellia and D. ananassae 


Rya-R44F 


3 


D. sechellia 


Nup98 


13 


D. melanogaster 


Nup214 


10 


D. erecta 


Nup50 


13 


D. ananassae 


CG6540 


13 


D. ananassae 



Note. — Column 2 listed the number of interactions of each gene within Nups 
genes network. 



Finally, it is worth to mention that both maternal and 
early zygotic genes evolved at similar rates in D. yakuba 
and D. erecta (fig. 1). However, these results rest on the 
assumption that the timing of expression of embryonic 
genes is conserved across the entire phylogeny of the 
D. melanogaster group. As a matter of fact, studies of 
genes with sex-biased expression in D. melanogaster and 



D. ananassae have shown that about one-third of the 
genes have either gained or lost sex-biased expression in 
one species (Grath et al. 2009). These changes in the pat- 
terns of gene expression across two distantly related species 
have likely influenced the evolution of the so-called sex- 
biased genes. 

PS in the Nuclear Pore Gene Network and Implications 
for Speciation 

Despite the pervasive role played by negative selection affect- 
ing genes expressed during embryogenesis, we detected cases 
of PSG that are involved in the three embryonic stages studied 
and in the six species of the D. melanogaster group (supple- 
mentary table S3, Supplementary Material online). Moreover, 
the incidence of PS is similar, even greater in some cases, than 
in the sets of genes expressed in postembryonic and adult flies 
(table 3). Even though we found PSG and overrepresentation 
of some GO terms among fast evolving genes in all stages of 
embryogenesis (table 2), only the maternal fraction contained 
a network of fast evolving genes that consist of a cluster of 
nuclear pore (Nups) genes. Nups encode components of the 
nuclear pore complex that form the channels that allow the 
transport of proteins and RNAs from the nucleus to the cyto- 
plasm and vice versa (Allen et al. 2000; Devos et al. 2006; Tran 
and Wente 2006). Interestingly, comparative genomics studies 
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Fig. 4. — Comparisons of the median 6N/65 values between Mvps and the rest of embryonic genes for the six species of the D. melanogaster group. Red 
asterisks show a significant difference in d/V/d5 between A/i/ps and rest of embryonic genes in a Kolmogorov-Smirnov test, P< 0.05. Black circles indicate 
ancestral nodes where a comparison between inferred sequences of the respective most recent common ancestors and derived lineages was performed. Red 
line shows an acceleration of nonsynonymous substitution rate of Nups in D. sechellia linage in comparison with its last common ancestor with D. simulans 
(Wilcoxon Matched Pairs test: N= 15, Z=2.7, P< 0.01). The rest of comparisons are not significant. 




Fig. 5. — Relative expression levels of the Nups network genes. All expression data were taken from Graveley et al. (201 1). Tau values, r, expressed 
temporal specificity of Nups genes across development. 
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indicate that a core of interacting proteins of the nuclear pore 
have been preserved for at least 1 .5 billion years, their associ- 
ation being at least as ancient as the last eukaryotic common 
ancestor (Bapteste et al. 2005; Neumann et al. 2010). Despite 
such ancient conservation, and in agreement with many 
reports, we found that Nups are fast evolving genes in 
Drosophila (Bapteste et al. 2005; Presgraves and Stephan 
2007; Tang and Presgraves 2009; Clark and Aquadro 2010). 
In any case, rapid evolution of Nups is at odds with the expec- 
tation that proteins involved in so relevant cellular mechanisms 
ought to be highly constrained and under negative selection. 
Several hypotheses have been proposed to explain Nups rapid 
evolution. On the one hand, Presgraves (2007) suggested 
that accelerated evolution of Nups may be related to nuclear 
transport-related segregation distortion. On the other hand, 
Sawamura et al. (2010) argued that rapid evolution of Nups 
may reflect genetic conflicts involving the nuclear entry of 
retroviruses and retrotransposons. However, it is difficult to 
reconcile these hypotheses with the novel evidence: The hall- 
mark of PS in Nups is only evident in D. sechellia (present 
paper) and D. mauritiana (Nolte et al. 2013) which are part 
of the D. simulans clade, a triad of very recently diverged 
species (Garrigan et al. 2012) (fig. 4). By contrast, the conflict 
over nuclear transport-related segregation distortion is 
thought to be an ancient genetic conflict even predating the 
D. melanogaster and D. simulans split (Presgraves and Stephan 
2007). Likewise, it is difficult to envisage how nuclear entry of 
retroviruses and retrotransposons would impose a lineage- 
specific acceleration in Nups only in the recently diverged 
species. Instead, our proposal is that such lineage-specific ac- 
celeration occurring in a short evolutionary timescale is likely a 
molecular signature of a reproductive isolation-related process 
framed in the context of early development as it is suggested 
by their highest expression level during embryogenesis (fig. 5). 
In this sense, even if we cannot rule out the possibility that 
Nups acceleration is a consequence of their function during 
nonembryonic stages of development, the fact that Nups 
genes present highest expression level during embryogenesis 
is a reliable indicator of an important embryonic function. 
Thus, even in the face of pleiotropy our results point out 
that Nups rapid evolution is related to their expression and 
role during early development. This explanation has, in addi- 
tion to the pattern of Nups evolution presented here, empirical 
and theoretical support. First, it has been shown that many 
Nups are involved in hybrid incompatibilities between pairs of 
species of the D. melanogaster subgroup, a feature that places 
Nups in the selected group of genes involved in early stages of 
speciation or "speciation genes" (Tang and Presgraves 2009; 
Sawamura et al. 2010). As proposed in the Dobzhansky- 
Muller model of evolution of postzyotic barriers to gene 
flow, independent adaptive fixations in diverging populations 
can lead to hybrid incompatibilities between interacting genes 
due to negative epistasis (reviewed in Coyne and Orr 2004). 
Indeed, coevolution among Nups (Clark and Aquadro 2010) 



can exacerbate the establishment of hybrid incompatibilities as 
a consequence of a "contagious" effect of interacting genes, 
because each substitution in a Nup would trigger a revolu- 
tionary episode of change among other components of the 
gene network and extend the occurrence of negative epista- 
sis. Second, Nups are involved in transcriptional regulation of a 
key reproductive trait in early development (Mendjan et al. 
2006; Mason and Goldfarb 2009), because the nuclear 
pore complex provides docking sites for chromatin (Kohler 
and Hurt 2010) and interacts with the X chromosome as 
part of the dosage compensation complex (Mendjan et al. 
2006). This is particularly interesting because dosage compen- 
sation in Drosophila takes place in early development as a key 
step of male sex determination (Bernstein and Cline 1994; 
Manu et al. 2013), and several studies have suggested that 
F1 hybrid lethality in crosses between D. melanogaster and 
D. simulans is due to dosage compensation failure (Orr 1 989). 
Thus, several features of this scenario lead us to propose that 
postzygotic isolation between species of the melanogaster 
subgroup might be partly the result of an improper dosage 
compensation caused by Nups functional divergence as it also 
happens for other components of the dosage compensation 
complex (Rodriguez et al. 2007; Sawamura 2012). If this 
hypothesis is correct, our study may help to understand 
how rapid evolving genes involved in the determination of a 
key reproductive trait in early development affect species 
divergence. 

Supplementary Material 

Supplementary tables S1-S3 are available at Genome Biology 
and Evolution online (http://www.gbe.oxfordjournals.org/). 
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